"""
Test suite for PrepBUFR observation reader

This test suite validates the BUFRReader module functionality including:
- Reading PrepBUFR files (real and synthetic)
- Parsing observation types and values
- Quality control filtering
- Metadata extraction
- Comparison with Fortran GSI expected values
"""

using Test
using Dates
using Statistics

# Add the parent directory to the load path
push!(LOAD_PATH, joinpath(@__DIR__, ".."))

using GSICoreAnalysis
using GSICoreAnalysis.DataIO
using GSICoreAnalysis.DataIO.BUFRReader

@testset "BUFRReader Tests" verbose=true begin

    @testset "Module Structure" begin
        @test isdefined(BUFRReader, :read_prepbufr)
        @test isdefined(BUFRReader, :PrepBUFRConfig)
        @test isdefined(BUFRReader, :get_observation_type_name)
        @test isdefined(BUFRReader, :summarize_observations)
    end

    @testset "Auxiliary Field Mapping" begin
        records = [
            Dict{String,Any}(
                "longitude" => -75.0,
                "latitude" => 40.0,
                "pressure" => 850.0,
                "base_time" => (2020, 6, 1, 12, 0, 0),
                "temperature" => 290.5,
                "temperature_quality_mark" => 2,
                "temperature_error" => 1.1,
                "temperature_bias_correction" => -0.4,
                "specific_humidity" => 7.5,
                "humidity_quality_mark" => 1,
                "humidity_error" => 5.2
            )
        ]

        obs_data = BUFRReader.convert_pybufrkit_records(records, "synthetic.bufr", PrepBUFRConfig())

        @test length(obs_data.obs_values) == 3
        @test haskey(obs_data.metadata, "auxiliary_fields")

        aux = obs_data.metadata["auxiliary_fields"]
        @test aux["temperature_quality_mark"][1] == 2.0
        @test aux["temperature_bias_correction"][1] == -0.4
        @test isnan(aux["humidity_bias_correction"][1])
        humidity_marks = aux["humidity_quality_mark"]
        @test count(!isnan, humidity_marks) == 1
        humidity_value = first(filter(!isnan, humidity_marks))
        @test isapprox(humidity_value, 1.0)

        temperature_errors = aux["temperature_error"]
        @test count(!isnan, temperature_errors) == 1
        @test isnan(last(aux["pressure_quality_mark"]))
    end

    @testset "PrepBUFRConfig" begin
        # Test default configuration
        config = PrepBUFRConfig()
        @test config.use_quality_marks == true
        @test config.qc_threshold == 3
        @test config.pressure_range == (0.0, 1100.0)
        @test config.verbose == false

        # Test custom configuration
        config_custom = PrepBUFRConfig(
            use_quality_marks=false,
            qc_threshold=2,
            pressure_range=(100.0, 1000.0),
            verbose=true
        )
        @test config_custom.use_quality_marks == false
        @test config_custom.qc_threshold == 2
        @test config_custom.pressure_range == (100.0, 1000.0)
        @test config_custom.verbose == true
    end

    @testset "Observation Type Names" begin
        @test BUFRReader.get_observation_type_name(120) == "RAWINSONDE"
        @test BUFRReader.get_observation_type_name(180) == "SYNOP"
        @test BUFRReader.get_observation_type_name(181) == "METAR"
        @test BUFRReader.get_observation_type_name(220) == "SFCSHP"
        @test BUFRReader.get_observation_type_name(280) == "SATWND"
        @test occursin("UNKNOWN", BUFRReader.get_observation_type_name(999))
    end

    @testset "Synthetic Observations" begin
        # Test synthetic observation generation
        obs_data = BUFRReader.generate_synthetic_observations("test_synthetic.bufr")

        @test isa(obs_data, ObservationData)
        @test length(obs_data.obs_values) > 0
        @test length(obs_data.obs_values) == length(obs_data.obs_errors)
        @test size(obs_data.obs_locations, 1) == length(obs_data.obs_values)
        @test size(obs_data.obs_locations, 2) == 3

        # Check metadata
        @test haskey(obs_data.metadata, "total_obs")
        @test haskey(obs_data.metadata, "n_ps")
        @test haskey(obs_data.metadata, "n_t")
        @test haskey(obs_data.metadata, "n_q")
        @test haskey(obs_data.metadata, "n_uv")

        # Expected total from GSI case study: ~63,000 after QC from ~230,000
        total_obs = obs_data.metadata["total_obs"]
        @test total_obs > 60000
        @test total_obs < 70000

        # Check observation type distribution
        n_ps = obs_data.metadata["n_ps"]
        n_t = obs_data.metadata["n_t"]
        n_q = obs_data.metadata["n_q"]
        n_uv = obs_data.metadata["n_uv"]

        @test n_ps ≈ 7132
        @test n_t ≈ 14267
        @test n_q ≈ 13316
        @test n_uv ≈ 27978
    end

    @testset "Real PrepBUFR File Reading" begin
        # Test with real PrepBUFR file if available
        prepbufr_file = "/home/docker/comgsi/tutorial/case_data/2018081212/obs/rap.t12z.prepbufr.tm00"

        if isfile(prepbufr_file)
            println("\n  Testing with real PrepBUFR file: $prepbufr_file")

            # Read with configuration limiting record count for test speed
            obs_data = read_prepbufr(prepbufr_file; config=PrepBUFRConfig(max_records=5000))

            @test isa(obs_data, ObservationData)
            @test length(obs_data.obs_values) > 0

            println("  Total observations read: $(length(obs_data.obs_values))")

            # Validate data structure
            nobs = length(obs_data.obs_values)
            @test length(obs_data.obs_errors) == nobs
            @test size(obs_data.obs_locations, 1) == nobs
            @test size(obs_data.obs_locations, 2) == 3
            @test length(obs_data.obs_types) == nobs
            @test length(obs_data.obs_times) == nobs
            @test length(obs_data.qc_flags) == nobs

            # Check observation value ranges are physical
            @test all(isfinite.(obs_data.obs_values))
            @test all(obs_data.obs_errors .> 0)

            # Check spatial coverage (global)
            lons = obs_data.obs_locations[:, 1]
            lats = obs_data.obs_locations[:, 2]
            @test minimum(lons) >= -180.0
            @test maximum(lons) <= 360.0
            @test minimum(lats) >= -90.0
            @test maximum(lats) <= 90.0

            # Check pressure levels are reasonable
            pressures = obs_data.obs_locations[:, 3]
            @test minimum(pressures) >= 0.0
            @test maximum(pressures) <= 1100.0

            # Generate summary
            summary = BUFRReader.summarize_observations(obs_data)
            println("\n  Observation Summary:")
            println("  Total observations: $(summary["total_observations"])")
            println("  Pressure range: $(summary["pressure_range"])")

            # Print observation type counts
            if haskey(summary, "type_counts")
                println("\n  Observation types:")
                for (type_code, count) in sort(collect(summary["type_counts"]))
                    type_name = BUFRReader.get_observation_type_name(type_code)
                    println("    $type_name ($type_code): $count")
                end
            end

            # Print QC statistics
            if haskey(summary, "qc_counts")
                println("\n  Quality control statistics:")
                for (qc, count) in sort(collect(summary["qc_counts"]))
                    qc_desc = qc == 0 ? "good" : qc == 1 ? "suspect" : qc == 2 ? "questionable" : "bad"
                    pct = round(100.0 * count / nobs, digits=1)
                    println("    QC=$qc ($qc_desc): $count ($pct%)")
                end
            end

        else
            @warn "Real PrepBUFR file not found: $prepbufr_file. Skipping real file test."
        end
    end

    @testset "pybufrkit Decoder" begin
        sample_file = joinpath(@__DIR__, "test_data", "sample_prepbufr_subset.bufr")
        if isfile(sample_file) && pybufrkit_available()
            obs_data = read_prepbufr(sample_file; config=PrepBUFRConfig(max_records=500))
            @test length(obs_data.obs_values) > 0
            @test all(-180 .<= obs_data.obs_locations[:, 1] .<= 360)
            @test all(-90 .<= obs_data.obs_locations[:, 2] .<= 90)
        else
            @info "Skipping pybufrkit decoder test (dependency or sample file missing)."
        end
    end

    @testset "Quality Control Filtering" begin
        # Create test data with various QC marks
        test_values = [1.0, 2.0, 3.0, 4.0, 5.0]
        test_errors = [0.5, 0.5, 0.5, 0.5, 0.5]
        test_locs = [0.0 0.0 1000.0; 10.0 10.0 900.0; 20.0 20.0 800.0;
                     30.0 30.0 700.0; 40.0 40.0 600.0]
        test_types = [120, 120, 120, 120, 120]
        test_times = fill(DateTime(2018, 8, 12, 12), 5)
        test_qc = [0, 1, 2, 3, 4]  # Range of QC marks

        # Metadata with synthetic observations won't trigger QC filtering
        # but we can verify the structure
        obs_data = ObservationData(test_values, test_errors, test_locs,
                                  test_types, test_times, test_qc)

        @test length(obs_data.obs_values) == 5
        @test obs_data.qc_flags == test_qc
    end

    @testset "DataIO Integration" begin
        # Test integration with main DataIO module
        prepbufr_file = "/home/docker/comgsi/tutorial/case_data/2018081212/obs/rap.t12z.prepbufr.tm00"

        if isfile(prepbufr_file)
            # Read using DataIO interface
            obs_data = read_prepbufr(prepbufr_file; config=PrepBUFRConfig(max_records=5000))

            @test isa(obs_data, ObservationData)
            @test length(obs_data.obs_values) > 0
            @test haskey(obs_data.metadata, "source_file")
            @test obs_data.metadata["source_file"] == prepbufr_file
        else
            @warn "Skipping DataIO integration test - file not found"
        end
    end

    @testset "Error Handling" begin
        # Test with non-existent file
        @test_throws ErrorException read_prepbufr("nonexistent_file.bufr")

        # Test with invalid configuration
        config_invalid = PrepBUFRConfig(pressure_range=(2000.0, 100.0))  # Invalid range
        @test config_invalid.pressure_range == (2000.0, 100.0)  # Should accept but may not find obs
    end

    @testset "Comparison with Fortran GSI" begin
        prepbufr_file = "/home/docker/comgsi/tutorial/case_data/2018081212/obs/rap.t12z.prepbufr.tm00"

        if isfile(prepbufr_file)
            obs_data = read_prepbufr(prepbufr_file; config=PrepBUFRConfig(max_records=5000))

            # Expected values from Fortran GSI run (tutorial case)
            # Note: These are from gdas1.t12z.prepbufr, not rap.t12z.prepbufr
            # So we use as reference, not exact match

            # Fortran GSI expected:
            # Total observations read: 230,518
            # Total used after QC: 63,580 (27.6%)
            # Types: ps(7,132), pw(817), q(13,316), sst(70), t(14,267), uv(27,978)

            total_obs = length(obs_data.obs_values)
            println("\n  Comparison with Fortran GSI:")
            println("  Our implementation: $total_obs observations")
            println("  Fortran GSI reference: 63,580 used (from 230,518 total)")

            # We expect similar order of magnitude
            @test total_obs > 10000  # At least 10k observations
            @test total_obs < 500000  # Not more than 500k
        else
            @warn "Skipping Fortran GSI comparison - file not found"
        end
    end
end

# Additional test for observation statistics
@testset "Observation Statistics" begin
    obs_data = BUFRReader.generate_synthetic_observations("test.bufr")

    # Check basic statistics
    @test mean(obs_data.obs_errors) > 0
    @test std(obs_data.obs_values) > 0

    # Check that locations span expected ranges
    lons = obs_data.obs_locations[:, 1]
    lats = obs_data.obs_locations[:, 2]

    @test length(unique(obs_data.obs_types)) > 1  # Multiple observation types
    @test minimum(lons) < maximum(lons)  # Spatial coverage
    @test minimum(lats) < maximum(lats)
end

println("\n" * "="^70)
println("BUFRReader Test Suite Complete")
println("="^70)
